import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import sklearn as skl
import scipy as sci
import sympy as sp
#**note** this next bit is to hide warning due to upcoming changes to libaries which are not critical for the course
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

madagascar_control = pd.read_excel('2024_madagasscar_project_data.xlsx')

#Radius of Earth
r = 6378000

madagascar_mercator = madagascar_control
madagascar_mercator

#converting lon and lat from degree to radian
madagascar_mercator['Lon_rad_E'] = np.deg2rad(madagascar_control['Long_deg_E'])
madagascar_mercator['Lat_rad_N'] = np.deg2rad(madagascar_control['Lat_deg_N'])
madagascar_mercator

#Convert control point geographical coordinates to Mercator Projection grid coordinates.
madagascar_mercator['Mercator_Easting_Latitude'] = r * madagascar_control['Lon_rad_E']
madagascar_mercator['Mercator_Northing_Longitude'] = r * np.log(np.tan(np.pi/4 + madagascar_control['Lat_rad_N']/2))

madagascar_control_0_2=madagascar_mercator[['Mercator_Easting_Latitude','Mercator_Northing_Longitude']].iloc[[0,2]].copy()
madagascar_mercator

## negative pixel y values
madagascar_mercator['negative_y'] = -madagascar_mercator['pixel_y']

madagascar_control_0_2_pix =madagascar_control[['pixel_x','negative_y']].iloc[[0,2]].copy()

col1= (pd.concat([madagascar_control_0_2_pix['pixel_x'], madagascar_control_0_2_pix['negative_y']]).sort_index(kind='merge').to_numpy())
print(col1)

col2 =(pd.concat([-madagascar_control_0_2_pix['negative_y'], madagascar_control_0_2_pix['pixel_x']]).sort_index(kind='merge').to_numpy())
print(col2)

col_tx = np.empty((len(col1),),float)
col_tx[::2] = 1
col_tx[1::2] = 0
print(col_tx)

col_ty = np.empty((len(col1),),float)
col_ty[::2] = 0
col_ty[1::2] = 1
print(col_ty)

[ 419 -163  586 -340]
[163 419 340 586]
[1. 0. 1. 0.]
[0. 1. 0. 1.]

a_0_2=np.column_stack((col1,col2,col_tx,col_ty))
print(a_0_2)

[[ 419.  163.    1.    0.]
 [-163.  419.    0.    1.]
 [ 586.  340.    1.    0.]
 [-340.  586.    0.    1.]]

inv_a_0_2 = np.linalg.inv(a_0_2)

print(inv_a_0_2)

[[-0.00282009  0.00298896  0.00282009 -0.00298896]
 [-0.00298896 -0.00282009  0.00298896  0.00282009]
 [ 2.66881691 -0.79269817 -1.66881691  0.79269817]
 [ 0.79269817  2.66881691 -0.79269817 -1.66881691]]

madagascar_control_0_2_mer =madagascar_control[['Mercator_Easting_Latitude','Mercator_Northing_Longitude']].iloc[[0,2]].copy()

mer_e=madagascar_control_0_2_mer['Mercator_Easting_Latitude']
mer_n=madagascar_control_0_2_mer['Mercator_Northing_Longitude']
l_vect = (pd.concat([mer_e,mer_n]).sort_index(kind='merge').to_numpy())
l_vect

array([ 3.33951299e+06, -7.08100245e-10,  6.67902598e+06, -3.50347459e+06])

x_array = np.dot(inv_a_0_2, l_vect)
print(x_array)

[ 1.98894537e+04  1.01549243e+02 -5.01072065e+06  3.19943183e+06]

a,b,t_x,t_y=x_array

print(a,b,t_x,t_y)

19889.45374098237 101.54924306154862 -5010720.653324694 3199431.8269373374

madagascar_mercator['mercator_x_2c']=((a*madagascar_mercator['pixel_x'
                                ]-b*madagascar_mercator['pixel_y']+1*t_x+0*t_y))
madagascar_mercator['mercator_y_2c']=((a*madagascar_mercator['pixel_y'
                                ]-b*madagascar_mercator['pixel_x']+0*t_x+1*t_y)*(-1))
madagascar_mercator.to_excel('madagascar_mercator_transformed.xlsx', index=False)
madagascar_mercator

# Plot the scatterplot with adjusted transparency and point size
plot = sns.scatterplot(
    data=madagascar_mercator, 
    x='mercator_x_2c', 
    y='mercator_y_2c', 
    hue='crtl_pt_ID', 
    palette='coolwarm',  # A harmonious color palette for better contrast
    s=120,               # Increase point size for better visibility
    alpha=0.8            # Adjust transparency to improve overlapping point visibility
)

# Move the legend to the upper left
sns.move_legend(plot, 'upper left', bbox_to_anchor=(1, 1))

# Set title, axis labels, and font size
plot.set_title("Madagascar Points in Mercator Projection", fontsize=18, fontweight='bold')
plot.set_xlabel("Mercator X", fontsize=14)
plot.set_ylabel("Mercator Y", fontsize=14)

# Adjust tick parameters for readability
plot.tick_params(axis='both', which='major', labelsize=12)

# Set the aspect ratio to 'equal' for accurate scaling
plot.set_aspect('equal')

# Display the plot with tight layout to avoid label clipping
plt.tight_layout()
plt.show()

pixel_x_co = madagascar_control['pixel_x'].to_numpy()
pixel_y_co = madagascar_control['negative_y'].to_numpy()

pixel_x_co

array([419, 586, 586, 419, 525, 532, 528, 528, 517, 513, 494, 497, 499,
       510, 517, 525])

pixel_y_co

array([-163, -163, -340, -340, -233, -252, -253, -259, -294, -309, -301,
       -272, -256, -254, -242, -233])

east_transformed = a * pixel_x_co - b * pixel_y_co + t_x
north_transformed = b * pixel_x_co + a *pixel_y_co + t_y

east_transformed
north_transformed

array([       0.        ,    16958.72359128, -3503474.5885626 ,
       -3520433.31215388, -1381497.54210424, -1758686.31848148,
       -1778981.9691947 , -1898318.6916406 , -2595566.61424866,
       -2894314.61733564, -2737128.42302595, -2160029.61680828,
       -1841595.25846644, -1800699.30931079, -1561315.01971758,
       -1381497.54210424])

madagascar_transformed = madagascar_mercator

madagascar_transformed

madagascar_mercator = madagascar_transformed
madagascar_transformed['East_Transformed'] = east_transformed
madagascar_transformed['North_Transformed'] = north_transformed
madagascar_transformed1=madagascar_transformed
madagascar_transformed

# Set seaborn style for a clean background
sns.set(style="whitegrid", palette="muted")  # A cleaner background with grid lines

# Create the plot
plt.figure(figsize=(10, 8))

# Plot the scatterplot with adjusted transparency and point size
plot = sns.scatterplot(
    data=madagascar_transformed, 
    x='East_Transformed', 
    y='North_Transformed', 
    hue='crtl_pt_ID', 
    palette='coolwarm',  # A harmonious color palette for better contrast
    s=120,               # Increase point size for better visibility
    alpha=0.8            # Adjust transparency to improve overlapping point visibility
)
# Set title, axis labels, and font size
plot.set_title("Madagascar Transformation", fontsize=18, fontweight='bold')
plot.set_xlabel("East Transformed (Longitude)", fontsize=14)
plot.set_ylabel("North Transformed (Latitude)", fontsize=14)

# Adjust tick parameters for readability
plot.tick_params(axis='both', which='major', labelsize=12)

# Set the aspect ratio to 'equal' for accurate scaling
plot.set_aspect('equal')

# Display the plot with tight layout to avoid label clipping
plt.tight_layout()
plt.show()

geo_latitude_rad = np.pi / 2 - 2 * np.arctan(np.exp(-madagascar_transformed['North_Transformed'] / r))
geo_longitude_rad = madagascar_mercator['East_Transformed'] / r 
## formula Of page 44

madagascar_transformed['radians_longitude'] = geo_longitude_rad 
madagascar_transformed['radians_latitude'] = geo_latitude_rad

madagascar_transformed['geo_longitude'] = geo_longitude_rad * (180 / np.pi)  
madagascar_transformed['geo_latitude'] = geo_latitude_rad * (180 / np.pi) 
madagascar_transformed

# Plot the scatterplot with adjusted transparency and point size
plt.figure(figsize=(10, 8))
plot = sns.scatterplot(
    data=madagascar_transformed, 
    x='geo_longitude', 
    y='geo_latitude', 
    hue='crtl_pt_ID', 
    palette='coolwarm',  # A harmonious color palette for better contrast
    s=120,               # Increase point size for better visibility
    alpha=0.8            # Adjust transparency to improve overlapping point visibility
)

# Set title, axis labels, and font size
plot.set_title("Madagascar Location", fontsize=18, fontweight='bold')
plot.set_xlabel("Longitude (geo_longitude)", fontsize=14)
plot.set_ylabel("Latitude (geo_latitude)", fontsize=14)

# Adjust tick parameters for readability
plot.tick_params(axis='both', which='major', labelsize=12)

# Set the aspect ratio to 'equal' for accurate scaling
plot.set_aspect('equal')

# Display the plot with tight layout to avoid label clipping
plt.tight_layout()
plt.show()

# define the latitude of origin (equator) and two standard parallels in degrees.
phi_0 = 0
phi_1 = -10
phi_2 = -25

# Define the latitude of origin (equator) and two standard parallels in degrees.
# Convert degrees to radians for mathematical calculations.
phi_0rad = np.deg2rad(phi_0)
phi_1rad = np.deg2rad(phi_1)
phi_2rad = np.deg2rad(phi_2)

# Calculate the LCC projection scale factor n.
# n determines the convergence rate of meridians based on the two standard parallels.
lcc_n = np.log(np.cos(phi_1rad) / np.cos(phi_2rad)) / np.log( np.tan(np.pi / 4 + phi_2rad/ 2) / np.tan(np.pi / 4 + phi_1rad / 2))
lcc_n

np.float64(-0.301570626498758)

lcc_f = (np.cos(phi_1rad) * (np.tan(np.pi / 4 + phi_1rad / 2)**lcc_n))/lcc_n
lcc_f

np.float64(-3.443007923614173)

# ompute the projection radius rho_0 at the latitude of origin (φ_0).
# This is a key parameter for transforming geographic to projected coordinates.
rho_0=r * lcc_f/np.tan(np.pi/4+phi_0rad/2)**lcc_n
rho_0

np.float64(-21959504.536811195)

# Convert latitude and longitude (in radians) to LCC projection coordinates:
# Radial distance: Derived from the scaling coefficient and input latitude.
# Longitude difference: Scaled using the factor n.

madagascar_transformed['2C_LCC_latitude'] = r * lcc_f/np.tan(np.pi/4+madagascar_transformed['radians_latitude']/2)**lcc_n

madagascar_transformed['2C_LCC_longitude'] = lcc_n * (madagascar_transformed['radians_longitude']-phi_0rad)

madagascar_lambert_conical_conformal = madagascar_transformed

# 2C_LCC_Easting: Compute the X coordinate using the sine of the transformed longitude.
# 2C_LCC_Northin: Compute the Y coordinate by adjusting the radial distance with the cosine of the transformed longitude.

madagascar_lambert_conical_conformal['2C_LCC_Easting']= np.sin(madagascar_lambert_conical_conformal['2C_LCC_longitude']) * madagascar_lambert_conical_conformal['2C_LCC_latitude']
madagascar_lambert_conical_conformal['2C_LCC_Northing']= rho_0 - madagascar_lambert_conical_conformal['2C_LCC_latitude'] * np.cos(madagascar_lambert_conical_conformal['2C_LCC_longitude'])
madagascar_lambert_conical_conformal

# Plot the scatterplot with adjusted transparency and point size
plt.figure(figsize=(10, 8))
plot = sns.scatterplot(
    data=madagascar_lambert_conical_conformal, 
    x='2C_LCC_Easting', 
    y='2C_LCC_Northing', 
    hue='crtl_pt_ID', 
    palette='coolwarm',  # A harmonious color palette for better contrast
    s=180,               # Increase point size for better visibility
    alpha=0.8            # Adjust transparency to improve overlapping point visibility
)

# Set title, axis labels, and font size
plot.set_title("Madagascar LCC", fontsize=18, fontweight='bold')
plot.set_xlabel("Eastings (2C_LCC_Easting)", fontsize=14)
plot.set_ylabel("Northings (2C_LCC_Northing)", fontsize=14)

# Adjust tick parameters for readability
plot.tick_params(axis='both', which='major', labelsize=12)

# Set the aspect ratio to 'equal' for accurate scaling
plot.set_aspect('equal')

# Display the plot with tight layout to avoid label clipping
plt.tight_layout()
plt.show()

madagascar_control_0_4 =madagascar_control.iloc[[0,1,2,3]].copy()
madagascar_control_0_4

similary_4c_l_os = (pd.concat([madagascar_control_0_4['Mercator_Easting_Latitude'], madagascar_control_0_4['Mercator_Northing_Longitude']]).sort_index(kind='merge').to_numpy())
similary_4c_l_os_break =pd.concat([madagascar_control_0_4['Mercator_Easting_Latitude'], madagascar_control_0_4['Mercator_Northing_Longitude']])
print(similary_4c_l_os_break)

0    3.339513e+06
1    6.679026e+06
2    6.679026e+06
3    3.339513e+06
0   -7.081002e-10
1   -7.081002e-10
2   -3.503475e+06
3   -3.503475e+06
dtype: float64

similary_4c_l_os_break_sort=similary_4c_l_os_break.sort_index(kind='merge')
similary_4c_l_os_break_sort

0    3.339513e+06
0   -7.081002e-10
1    6.679026e+06
1   -7.081002e-10
2    6.679026e+06
2   -3.503475e+06
3    3.339513e+06
3   -3.503475e+06
dtype: float64

similary_4c_l_os_break_sort_numpy = similary_4c_l_os_break_sort.to_numpy()
similary_4c_l_os_break_sort_numpy

array([ 3.33951299e+06, -7.08100245e-10,  6.67902598e+06, -7.08100245e-10,
        6.67902598e+06, -3.50347459e+06,  3.33951299e+06, -3.50347459e+06])

similary_4c_col1= (pd.concat([madagascar_control_0_4['pixel_x'], madagascar_control_0_4['negative_y']]).sort_index(kind='merge').to_numpy())
print(similary_4c_col1)

similary_4c_col2 =(pd.concat([-madagascar_control_0_4['negative_y'], madagascar_control_0_4['pixel_x']]).sort_index(kind='merge').to_numpy())
print(similary_4c_col2)

similary_4c_col_tx = np.empty((len(similary_4c_col1),),float)
similary_4c_col_tx[::2] = 1
similary_4c_col_tx[1::2] = 0
print(similary_4c_col_tx)

similary_4c_col_ty = np.empty((len(similary_4c_col1),),float)
similary_4c_col_ty[::2] = 0
similary_4c_col_ty[1::2] = 1
print(similary_4c_col_ty)

[ 419 -163  586 -163  586 -340  419 -340]
[163 419 163 586 340 586 340 419]
[1. 0. 1. 0. 1. 0. 1. 0.]
[0. 1. 0. 1. 0. 1. 0. 1.]

similary_4c_stack=np.column_stack((similary_4c_col1,similary_4c_col2,similary_4c_col_tx,similary_4c_col_ty))
print(similary_4c_stack)

[[ 419.  163.    1.    0.]
 [-163.  419.    0.    1.]
 [ 586.  163.    1.    0.]
 [-163.  586.    0.    1.]
 [ 586.  340.    1.    0.]
 [-340.  586.    0.    1.]
 [ 419.  340.    1.    0.]
 [-340.  419.    0.    1.]]

similary_4c_v_vect = np.matmul(similary_4c_stack.T, similary_4c_l_os_break_sort_numpy)

similary_4c_ata=np.matmul(similary_4c_stack.T, similary_4c_stack)
similary_4c_ata

array([[ 1.322252e+06,  0.000000e+00,  2.010000e+03, -1.006000e+03],
       [ 0.000000e+00,  1.322252e+06,  1.006000e+03,  2.010000e+03],
       [ 2.010000e+03,  1.006000e+03,  4.000000e+00,  0.000000e+00],
       [-1.006000e+03,  2.010000e+03,  0.000000e+00,  4.000000e+00]])

similary_4c_q_matrix = np.linalg.inv(similary_4c_ata)
print(similary_4c_q_matrix)

[[ 1.68867574e-05  1.34343290e-21 -8.48559560e-03  4.24701949e-03]
 [-1.34343290e-21  1.68867574e-05 -4.24701949e-03 -8.48559560e-03]
 [-8.48559560e-03 -4.24701949e-03  5.58213719e+00  0.00000000e+00]
 [ 4.24701949e-03 -8.48559560e-03  0.00000000e+00  5.58213719e+00]]

similary_4c_x_vect = np.matmul(similary_4c_q_matrix, similary_4c_v_vect)
print(similary_4c_x_vect)

[ 1.98894537e+04 -1.45519152e-11 -4.98518102e+06  3.25046032e+06]

## multiplication between A and X (E= A*X-L)

similary_4c_lin_model=np.matmul(similary_4c_stack,similary_4c_x_vect)

#lin_model=lin_model[:,np.newaxis]
#similary_4c_lin_model=similary_4c_lin_model
similary_4c_lin_model

array([ 3348500.0987769 ,     8479.36179563,  6670038.87352094,
           8479.36179563,  6670038.87352094, -3511953.95035824,
        3348500.0987769 , -3511953.95035824])

similary_4c_ls_errors=(similary_4c_lin_model - similary_4c_l_os_break_sort_numpy)
similary_4c_ls_errors=similary_4c_ls_errors.reshape(4,2)

x_ls = similary_4c_ls_errors[:, 0]
y_ls = similary_4c_ls_errors[:, 1]

madagascar_control_0_4['Mercat_x_ls_error'] = x_ls
madagascar_control_0_4['Mercat_y_ls_error'] = y_ls
madagascar_control_0_4

a_4c,b_4c,t_x_4C,t_y_4C=similary_4c_x_vect

similary_4c_east_transformed = a_4c * pixel_x_co - b_4c * pixel_y_co + t_x_4C
similary_4c_north_transformed = b_4c * pixel_x_co + a_4c *pixel_y_co + t_y_4C

Similary_4C_madagascar_transformed = madagascar_lambert_conical_conformal

Similary_4C_madagascar_transformed['4C_East_Mercator_Transformed'] = similary_4c_east_transformed
Similary_4C_madagascar_transformed['4C_North_Mercator_Transformed'] = similary_4c_north_transformed
Similary_4C_madagascar_transformed

# Plot the scatterplot with adjusted transparency and point size
plt.figure(figsize=(10, 8))
plot = sns.scatterplot(
    data=Similary_4C_madagascar_transformed, 
    x='4C_East_Mercator_Transformed', 
    y='4C_North_Mercator_Transformed', 
    hue='crtl_pt_ID', 
    palette='coolwarm',  # A harmonious color palette for better contrast
    s=100,               # Increase point size for better visibility
    alpha=0.8            # Adjust transparency to improve overlapping point visibility
)
# Set title, axis labels, and font size
plot.set_title("Madagascar Points Displayed in Mercator Projection (4 control points)", fontsize=11, fontweight='bold')
plot.set_xlabel("East", fontsize=14)
plot.set_ylabel("North", fontsize=14)

# Adjust tick parameters for readability
plot.tick_params(axis='both', which='major', labelsize=20)

# Set the aspect ratio to 'equal' for accurate scaling
plot.set_aspect('equal')

# Display the plot with tight layout to avoid label clipping
plt.tight_layout()
plt.show()

## formula on page 44
similary_4c_geo_latitude_rad = np.pi / 2 - 2 * np.arctan(np.exp(-Similary_4C_madagascar_transformed['4C_North_Mercator_Transformed'] / r))
similary_4c_geo_longitude_rad =Similary_4C_madagascar_transformed['4C_East_Mercator_Transformed'] / r

similary_4c_geo_longitude = similary_4c_geo_longitude_rad * (180 / np.pi)  
similary_4c_geo_latitude = similary_4c_geo_latitude_rad * (180 / np.pi)

madagascar_transformed['4C_geo_longitude'] = similary_4c_geo_longitude
madagascar_transformed['4C_geo_latitude'] = similary_4c_geo_latitude
madagascar_transformed

Similary_4C_madagascar_transformed['4C_geo_latitude'] = r * lcc_f/np.tan(np.pi/4+similary_4c_geo_latitude_rad /2)**lcc_n

Similary_4C_madagascar_transformed['4C_geo_longtitude'] = lcc_n * (similary_4c_geo_longitude_rad-phi_0rad)

Similary_4C_madagascar_transformed['4C_LCC_Easting']= np.sin(Similary_4C_madagascar_transformed['4C_geo_longtitude']) * Similary_4C_madagascar_transformed['4C_geo_latitude']
Similary_4C_madagascar_transformed['4C_LCC_Northing']= rho_0 - Similary_4C_madagascar_transformed['4C_geo_latitude'] * np.cos(Similary_4C_madagascar_transformed['4C_geo_longtitude'])
Similary_4C_madagascar_transformed

# Plot the scatterplot with adjusted transparency and point size
plt.figure(figsize=(10, 8))
plot = sns.scatterplot(
    data=Similary_4C_madagascar_transformed, 
    x='4C_LCC_Easting', 
    y='4C_LCC_Northing', 
    hue='crtl_pt_ID', 
    palette='coolwarm',  # A harmonious color palette for better contrast
    s=150,               # Increase point size for better visibility
    alpha=0.8            # Adjust transparency to improve overlapping point visibility
)

# Set title, axis labels, and font size
plot.set_title("Madagascar Points Displayed in Lambert Conformal Conic Projection (4 control points)", fontsize=10, fontweight='bold')
plot.set_xlabel("Eastings (4C_LCC_Easting)", fontsize=8)
plot.set_ylabel("Northings (4C_LCC_Northing)", fontsize=8)

# Adjust tick parameters for readability
plot.tick_params(axis='both', which='major', labelsize=8)

# Set the aspect ratio to 'equal' for accurate scaling
plot.set_aspect('equal')

# Display the plot with tight layout to avoid label clipping
plt.tight_layout()
plt.show()

Similary_4C_madagascar_transformed

#madagascar_control_0_4['Least Squares Fitting Mercator easting(x)']=similary_4c_lin_model[:,0]
#madagascar_control_0_4['Least Squares Fitting Mercator northing(y)']=similary_4c_lin_model[:,1]

Similary_4C_madagascar_transformed['Least Squares Fitting LCC residual'] =np.sqrt((Similary_4C_madagascar_transformed['4C_LCC_Easting'])**2 + (Similary_4C_madagascar_transformed['4C_LCC_Northing'])**2)
Similary_4C_madagascar_transformed

# Plot the scatterplot with adjusted transparency and point size
plt.figure(figsize=(10, 8))
plot = sns.scatterplot(
    data=madagascar_lambert_conical_conformal, 
    x='4C_LCC_Easting', 
    y='4C_LCC_Northing', 
    hue='crtl_pt_ID', 
    palette='coolwarm',  # A harmonious color palette for better contrast
    s=180,               # Increase point size for better visibility
    alpha=0.8            # Adjust transparency to improve overlapping point visibility
)

# Set title, axis labels, and font size
plot.set_title("Madagascar Lambert Conical Conformal (4 Control Points)", fontsize=18, fontweight='bold')
plot.set_xlabel("Eastings (4C_LCC_Easting)", fontsize=14)
plot.set_ylabel("Northings (4C_LCC_Northing)", fontsize=14)

# Adjust tick parameters for readability
plot.tick_params(axis='both', which='major', labelsize=12)

# Set the aspect ratio to 'equal' for accurate scaling
plot.set_aspect('equal')

# Display the plot with tight layout to avoid label clipping
plt.tight_layout()
plt.show()

madagascar_control_0_4

important_columns = ['crtl_pt_ID', 'Lat_deg_N','Long_deg_E', 'pixel_x','pixel_y', 'Lon_rad_E','Lat_rad_N', 'Mercator_Easting_Latitude','Mercator_Northing_Longitude', 'negative_y','2C_LCC_Easting', '2C_LCC_Northing','4C_East_Mercator_Transformed', '4C_North_Mercator_Transformed','4C_LCC_Easting', '4C_LCC_Northing']  

processed_4c_madagascar_transformed = Similary_4C_madagascar_transformed[important_columns]
processed_4c_madagascar_transformed

flight_coordinates= madagascar_transformed1.iloc[[8,11,15]].copy()
flight_coordinates

point_8_latitude = flight_coordinates['geo_latitude'].iloc[1]
point_8_longitude = flight_coordinates['geo_longitude'].iloc[1]

point_8_rad_latitude = np.radians(point_8_latitude)
point_8_rad_longitude = np.radians(point_8_longitude)

point_12_latitude = flight_coordinates['geo_latitude'].iloc[2] 
point_12_longitude = flight_coordinates['geo_longitude'].iloc[2]
point_12_rad_latitude = np.radians(point_12_latitude)
point_12_rad_longitude = np.radians(point_12_longitude)

point_12_rad_latitude = np.radians(flight_coordinates['geo_latitude'].iloc[2]) 
point_12_rad_longitude = np.radians(flight_coordinates['geo_longitude'].iloc[2])

cos_n = (np.sin(point_8_rad_latitude) * np.sin(point_12_rad_latitude) +
         np.cos(point_8_rad_latitude) * np.cos(point_12_rad_latitude) *
         np.cos(point_12_rad_longitude - point_8_rad_longitude))

n_radians = np.arccos(cos_n)
n_radians
n_degrees = np.rad2deg(n_radians)
_arc = np.array((n_radians, n_degrees))
round_arc= np.round(_arc,4)

earth_radius_km = 6378
distance_12_8_km = earth_radius_km * n_radians

point_5_latitude = flight_coordinates['geo_latitude'].iloc[0]
point_5_longitude = flight_coordinates['geo_longitude'].iloc[0]

point_5_rad_latitude = np.radians(flight_coordinates['geo_latitude'].iloc[0])
point_5_rad_longitude = np.radians(flight_coordinates['geo_longitude'].iloc[0])

cos_n = (np.sin(point_5_rad_latitude) * np.sin(point_8_rad_latitude) +
         np.cos(point_5_rad_latitude) * np.cos(point_8_rad_latitude) *
         np.cos(point_8_rad_longitude - point_5_rad_longitude))

n_radians = np.arccos(cos_n)
n_radians
n_degrees = np.rad2deg(n_radians)
_arc=np.array((n_radians,n_degrees))
round_arc= np.round(_arc,4)

distance_8_5_km = earth_radius_km * n_radians

total_distance=distance_12_8_km + distance_8_5_km
print(str(total_distance) + ' km')

1471.1765624218322 km

point_8_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[1]
point_8_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[1]
point_12_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[2]
point_12_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[2]
distance_12_8_mercator_2c = np.sqrt((point_12_mercator_longitude - point_8_mercator_longitude)**2 +
                                    (point_12_mercator_latitude - point_8_mercator_latitude)**2)
distance_12_8_mercator_km_2c = distance_12_8_mercator_2c / 1000
point_5_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[0]
point_5_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[0]
distance_8_5_mercator_2c = np.sqrt((point_8_mercator_longitude - point_5_mercator_longitude)**2 +
                                    (point_8_mercator_latitude - point_5_mercator_latitude)**2)
distance_8_5_mercator_km_2c = distance_8_5_mercator_2c / 1000
total_mercator_projection_distance_2c= distance_12_8_mercator_km_2c + distance_8_5_mercator_km_2c
print(str(total_mercator_projection_distance_2c) + ' km')

1547.8721224618705 km

# Extract coordinates for the points
point_8_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[1]
point_8_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[1]
point_12_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[2]
point_12_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[2]
point_5_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[0]
point_5_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[0]

# Calculate distance between C12 and C8
distance_12_8_mercator = np.sqrt((point_12_mercator_longitude - point_8_mercator_longitude)**2 +
                                 (point_12_mercator_latitude - point_8_mercator_latitude)**2)
distance_12_8_mercator_km = distance_12_8_mercator / 1000  # Convert to kilometers

# Calculate distance between C8 and C5
distance_8_5_mercator = np.sqrt((point_8_mercator_longitude - point_5_mercator_longitude)**2 +
                                (point_8_mercator_latitude - point_5_mercator_latitude)**2)
distance_8_5_mercator_km = distance_8_5_mercator / 1000  # Convert to kilometers

# Calculate total distance
total_mercator_projection_distance = distance_12_8_mercator_km + distance_8_5_mercator_km
print(str(total_mercator_projection_distance) + ' km')

# Plot the flight route
plt.figure(figsize=(8, 6))

# Plot the route from C12 to C8
plt.plot([point_12_mercator_longitude, point_8_mercator_longitude],
         [point_12_mercator_latitude, point_8_mercator_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8")

# Plot the route from C8 to C5
plt.plot([point_8_mercator_longitude, point_5_mercator_longitude],
         [point_8_mercator_latitude, point_5_mercator_latitude],
         color='red', marker='o', markersize=8, label="C8 to C5")

# Annotate the points
plt.text(point_12_mercator_longitude, point_12_mercator_latitude, 'C12', fontsize=12, ha='right')
plt.text(point_8_mercator_longitude, point_8_mercator_latitude, 'C8', fontsize=12, ha='right')
plt.text(point_5_mercator_longitude, point_5_mercator_latitude, 'C5', fontsize=12, ha='right')

# Set the title and labels
plt.title("Flight Route-Mercator distance(2 controls point)")
plt.xlabel("Longitude")
plt.ylabel("Latitude")

# Show the legend
plt.legend()

# Display the plot
plt.tight_layout()
plt.show()

1547.8721224618705 km

point_8_mercator_latitude = flight_coordinates['North_Transformed'].iloc[1]
point_8_mercator_longitude = flight_coordinates['East_Transformed'].iloc[1]

point_12_mercator_latitude = flight_coordinates['North_Transformed'].iloc[2]
point_12_mercator_longitude = flight_coordinates['East_Transformed'].iloc[2]

distance_12_8_mercator = np.sqrt((point_12_mercator_longitude - point_8_mercator_longitude)**2 +
                                    (point_12_mercator_latitude - point_8_mercator_latitude)**2)

distance_12_8_mercator_km = distance_12_8_mercator / 1000

point_5_mercator_latitude = flight_coordinates['North_Transformed'].iloc[0]
point_5_mercator_longitude = flight_coordinates['East_Transformed'].iloc[0]

distance_8_5_mercator = np.sqrt((point_8_mercator_longitude - point_5_mercator_longitude)**2 +
                                    (point_8_mercator_latitude - point_5_mercator_latitude)**2)

distance_8_5_mercator_km = distance_8_5_mercator / 1000

total_mercator_projection_distance=distance_12_8_mercator_km + distance_8_5_mercator_km
print(str(total_mercator_projection_distance) + ' km')

1546.2770660005715 km

# Extract coordinates for the points from the transformed flight coordinates
point_8_mercator_latitude = flight_coordinates['North_Transformed'].iloc[1]
point_8_mercator_longitude = flight_coordinates['East_Transformed'].iloc[1]
point_12_mercator_latitude = flight_coordinates['North_Transformed'].iloc[2]
point_12_mercator_longitude = flight_coordinates['East_Transformed'].iloc[2]
point_5_mercator_latitude = flight_coordinates['North_Transformed'].iloc[0]
point_5_mercator_longitude = flight_coordinates['East_Transformed'].iloc[0]

# Calculate distance between C12 and C8
distance_12_8_mercator = np.sqrt((point_12_mercator_longitude - point_8_mercator_longitude)**2 +
                                 (point_12_mercator_latitude - point_8_mercator_latitude)**2)
distance_12_8_mercator_km = distance_12_8_mercator / 1000  # Convert to kilometers

# Calculate distance between C8 and C5
distance_8_5_mercator = np.sqrt((point_8_mercator_longitude - point_5_mercator_longitude)**2 +
                                (point_8_mercator_latitude - point_5_mercator_latitude)**2)
distance_8_5_mercator_km = distance_8_5_mercator / 1000  # Convert to kilometers

# Calculate total distance
total_mercator_projection_distance = distance_12_8_mercator_km + distance_8_5_mercator_km
print(f'Total distance: {total_mercator_projection_distance} km')

# Plot the flight route
plt.figure(figsize=(8, 6))

# Plot the route from C12 to C8
plt.plot([point_12_mercator_longitude, point_8_mercator_longitude],
         [point_12_mercator_latitude, point_8_mercator_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8")

# Plot the route from C8 to C5
plt.plot([point_8_mercator_longitude, point_5_mercator_longitude],
         [point_8_mercator_latitude, point_5_mercator_latitude],
         color='red', marker='o', markersize=8, label="C8 to C5")

# Annotate the points on the plot
plt.text(point_12_mercator_longitude, point_12_mercator_latitude, 'C12', fontsize=12, ha='right')
plt.text(point_8_mercator_longitude, point_8_mercator_latitude, 'C8', fontsize=12, ha='right')
plt.text(point_5_mercator_longitude, point_5_mercator_latitude, 'C5', fontsize=12, ha='right')

# Set the title and labels for the plot
plt.title("Flight Route-Mercator distance(4 controls point)")
plt.xlabel("Longitude")
plt.ylabel("Latitude")

# Display the legend
plt.legend()

# Show the plot
plt.tight_layout()
plt.show()

Total distance: 1546.2770660005715 km

point_8_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[1]
point_8_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[1]

pt12_l_lat = flight_coordinates['2C_LCC_Northing'].iloc[2]
pt12_l_long = flight_coordinates['2C_LCC_Easting'].iloc[2]

dist_12_8_l = np.sqrt((pt12_l_long - point_8_lambert_longitude)**2 +
                                    (pt12_l_lat - point_8_lambert_latitude)**2)

dist_12_8_l_km = dist_12_8_l / 1000

print('the Lambert Conformal Conic projection distance between 12 and 8 is {} km'.format(np.round(dist_12_8_l_km,1)))

the Lambert Conformal Conic projection distance between 12 and 8 is 911.9 km

point_5_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[0]
point_5_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[0]

distance_8_5_lambert = np.sqrt((point_8_lambert_longitude - point_5_lambert_longitude)**2 +
                                    (point_8_lambert_latitude - point_5_lambert_latitude)**2)

distance_8_5_lambert_km = distance_8_5_lambert/1000

total_l_dist=dist_12_8_l_km + distance_8_5_lambert_km
print(str(total_l_dist) + ' km')

1460.612996865721 km

# Extract coordinates for the points in Lambert Conformal Conic (LCC) projection
point_8_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[1]
point_8_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[1]
pt12_l_lat = flight_coordinates['2C_LCC_Northing'].iloc[2]
pt12_l_long = flight_coordinates['2C_LCC_Easting'].iloc[2]
point_5_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[0]
point_5_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[0]

# Calculate the distance between C12 and C8 in Lambert Conformal Conic (LCC) coordinates
dist_12_8_l = np.sqrt((pt12_l_long - point_8_lambert_longitude)**2 +
                      (pt12_l_lat - point_8_lambert_latitude)**2)
dist_12_8_l_km = dist_12_8_l / 1000  # Convert to kilometers

# Calculate the distance between C8 and C5 in Lambert Conformal Conic (LCC) coordinates
distance_8_5_lambert = np.sqrt((point_8_lambert_longitude - point_5_lambert_longitude)**2 +
                                (point_8_lambert_latitude - point_5_lambert_latitude)**2)
distance_8_5_lambert_km = distance_8_5_lambert / 1000  # Convert to kilometers

# Calculate total distance from C12 to C8 to C5
total_lambert_projection_distance = dist_12_8_l_km + distance_8_5_lambert_km
print(f'Total distance from C12 to C8 to C5: {total_lambert_projection_distance} km')

# Plot the flight route
plt.figure(figsize=(8, 6))

# Plot the route from C12 to C8
plt.plot([pt12_l_long, point_8_lambert_longitude],
         [pt12_l_lat, point_8_lambert_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8")

# Plot the route from C8 to C5
plt.plot([point_8_lambert_longitude, point_5_lambert_longitude],
         [point_8_lambert_latitude, point_5_lambert_latitude],
         color='red', marker='o', markersize=8, label="C8 to C5")

# Annotate the points on the plot
plt.text(pt12_l_long, pt12_l_lat, 'C12', fontsize=12, ha='right')
plt.text(point_8_lambert_longitude, point_8_lambert_latitude, 'C8', fontsize=12, ha='right')
plt.text(point_5_lambert_longitude, point_5_lambert_latitude, 'C5', fontsize=12, ha='right')

# Set the title and labels for the plot
plt.title("Flight Route from C12 to C8 to C5 (Lambert Conformal Conic Projection)")
plt.xlabel("Easting (Longitude in Lambert Conformal Conic)")
plt.ylabel("Northing (Latitude in Lambert Conformal Conic)")

# Display the legend
plt.legend()

# Show the plot
plt.tight_layout()
plt.show()

Total distance from C12 to C8 to C5: 1460.612996865721 km

point_8_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[1]
point_8_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[1]

point_12_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[2]
point_12_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[2]

distance_12_8_4c_lambert = np.sqrt((point_12_4c_lambert_longitude - point_8_4c_lambert_longitude)**2 +
                                    (point_12_4c_lambert_latitude - point_8_4c_lambert_latitude)**2)

distance_12_8_4c_lambert_km = distance_12_8_4c_lambert / 1000

print('the Lambert Conformal Conic projection distance between 12 and 8 is {} km (4 control points)'.format(np.round(distance_12_8_4c_lambert_km,1)))

the Lambert Conformal Conic projection distance between 12 and 8 is 911.8 km (4 control points)

point_5_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[0]
point_5_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[0]

distance_8_5_4c_lambert = np.sqrt((point_8_4c_lambert_longitude - point_5_4c_lambert_longitude)**2 +
                                    (point_8_4c_lambert_latitude - point_5_4c_lambert_latitude)**2)

distance_8_5_4c_lambert_km = distance_8_5_4c_lambert/1000

total_lambert_projection_4c_distance=distance_12_8_4c_lambert_km + distance_8_5_4c_lambert_km
print(str(total_lambert_projection_4c_distance) + ' km')

1460.5443208825818 km

# Extract coordinates for the points in the 4C Lambert Conformal Conic (LCC) projection
point_8_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[1]
point_8_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[1]
point_12_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[2]
point_12_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[2]
point_5_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[0]
point_5_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[0]

# Calculate the distance between C12 and C8 in 4C Lambert Conformal Conic (LCC) coordinates
distance_12_8_4c_lambert = np.sqrt((point_12_4c_lambert_longitude - point_8_4c_lambert_longitude)**2 +
                                   (point_12_4c_lambert_latitude - point_8_4c_lambert_latitude)**2)
distance_12_8_4c_lambert_km = distance_12_8_4c_lambert / 1000  # Convert to kilometers

# Calculate the distance between C8 and C5 in 4C Lambert Conformal Conic (LCC) coordinates
distance_8_5_4c_lambert = np.sqrt((point_8_4c_lambert_longitude - point_5_4c_lambert_longitude)**2 +
                                  (point_8_4c_lambert_latitude - point_5_4c_lambert_latitude)**2)
distance_8_5_4c_lambert_km = distance_8_5_4c_lambert / 1000  # Convert to kilometers

# Calculate total distance from C12 to C8 to C5
total_4c_lambert_projection_distance = distance_12_8_4c_lambert_km + distance_8_5_4c_lambert_km
print(f'Total distance from C12 to C8 to C5: {total_4c_lambert_projection_distance} km')

# Plot the flight route using Lambert Conformal Conic (LCC) coordinates
plt.figure(figsize=(8, 6))

# Plot the route from C12 to C8
plt.plot([point_12_4c_lambert_longitude, point_8_4c_lambert_longitude],
         [point_12_4c_lambert_latitude, point_8_4c_lambert_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8")

# Plot the route from C8 to C5
plt.plot([point_8_4c_lambert_longitude, point_5_4c_lambert_longitude],
         [point_8_4c_lambert_latitude, point_5_4c_lambert_latitude],
         color='red', marker='o', markersize=8, label="C8 to C5")

# Annotate the points on the plot
plt.text(point_12_4c_lambert_longitude, point_12_4c_lambert_latitude, 'C12', fontsize=12, ha='right')
plt.text(point_8_4c_lambert_longitude, point_8_4c_lambert_latitude, 'C8', fontsize=12, ha='right')
plt.text(point_5_4c_lambert_longitude, point_5_4c_lambert_latitude, 'C5', fontsize=12, ha='right')

# Set the title and labels for the plot
plt.title("Flight Route from C12 to C8 to C5 (4C Lambert Conformal Conic Projection)")
plt.xlabel("Easting (Longitude in Lambert Conformal Conic)")
plt.ylabel("Northing (Latitude in Lambert Conformal Conic)")

# Display the legend
plt.legend()

# Show the plot
plt.tight_layout()
plt.show()

Total distance from C12 to C8 to C5: 1460.5443208825818 km

# 2C Mercator Projection (First Set of Coordinates)
point_8_2C_latitude = flight_coordinates['mercator_y_2c'].iloc[1]
point_8_2C_longitude = flight_coordinates['mercator_x_2c'].iloc[1]
point_12_2C_latitude = flight_coordinates['mercator_y_2c'].iloc[2]
point_12_2C_longitude = flight_coordinates['mercator_x_2c'].iloc[2]
point_5_2C_latitude = flight_coordinates['mercator_y_2c'].iloc[0]
point_5_2C_longitude = flight_coordinates['mercator_x_2c'].iloc[0]

# 4C Mercator Projection (Second Set of Coordinates)
point_8_4C_latitude = flight_coordinates['North_Transformed'].iloc[1]
point_8_4C_longitude = flight_coordinates['East_Transformed'].iloc[1]
point_12_4C_latitude = flight_coordinates['North_Transformed'].iloc[2]
point_12_4C_longitude = flight_coordinates['East_Transformed'].iloc[2]
point_5_4C_latitude = flight_coordinates['North_Transformed'].iloc[0]
point_5_4C_longitude = flight_coordinates['East_Transformed'].iloc[0]

# Calculate Distance for 2C Mercator Projection
distance_12_8_mercator_2C = np.sqrt((point_12_2C_longitude - point_8_2C_longitude)**2 +
                                    (point_12_2C_latitude - point_8_2C_latitude)**2)
distance_12_8_mercator_km_2C = distance_12_8_mercator_2C / 1000

distance_8_5_mercator_2C = np.sqrt((point_8_2C_longitude - point_5_2C_longitude)**2 +
                                    (point_8_2C_latitude - point_5_2C_latitude)**2)
distance_8_5_mercator_km_2C = distance_8_5_mercator_2C / 1000

total_mercator_projection_distance_2C = distance_12_8_mercator_km_2C + distance_8_5_mercator_km_2C
print(f"Total distance 2C Mercator: {total_mercator_projection_distance_2C} km")

# Calculate Distance for 4C Mercator Projection
distance_12_8_mercator_4C = np.sqrt((point_12_4C_longitude - point_8_4C_longitude)**2 +
                                    (point_12_4C_latitude - point_8_4C_latitude)**2)
distance_12_8_mercator_km_4C = distance_12_8_mercator_4C / 1000

distance_8_5_mercator_4C = np.sqrt((point_8_4C_longitude - point_5_4C_longitude)**2 +
                                    (point_8_4C_latitude - point_5_4C_latitude)**2)
distance_8_5_mercator_km_4C = distance_8_5_mercator_4C / 1000

total_mercator_projection_distance_4C = distance_12_8_mercator_km_4C + distance_8_5_mercator_km_4C
print(f"Total distance 4C Mercator: {total_mercator_projection_distance_4C} km")

# Plot the flight route for both projections
plt.figure(figsize=(8, 6))

# Plot the route from C12 to C8 for 2C Mercator
plt.plot([point_12_2C_longitude, point_8_2C_longitude],
         [point_12_2C_latitude, point_8_2C_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8 (2C)")

# Plot the route from C8 to C5 for 2C Mercator
plt.plot([point_8_2C_longitude, point_5_2C_longitude],
         [point_8_2C_latitude, point_5_2C_latitude],
         color='red', marker='o', markersize=8, label="C8 to C5 (2C)")

# Plot the route from C12 to C8 for 4C Mercator
plt.plot([point_12_4C_longitude, point_8_4C_longitude],
         [point_12_4C_latitude, point_8_4C_latitude],
         color='green', marker='o', markersize=8, label="C12 to C8 (4C)")

# Plot the route from C8 to C5 for 4C Mercator
plt.plot([point_8_4C_longitude, point_5_4C_longitude],
         [point_8_4C_latitude, point_5_4C_latitude],
         color='purple', marker='o', markersize=8, label="C8 to C5 (4C)")

# Annotate the points
plt.text(point_12_2C_longitude, point_12_2C_latitude, 'C12 (2C)', fontsize=12, ha='right')
plt.text(point_8_2C_longitude, point_8_2C_latitude, 'C8 (2C)', fontsize=12, ha='right')
plt.text(point_5_2C_longitude, point_5_2C_latitude, 'C5 (2C)', fontsize=12, ha='right')
plt.text(point_12_4C_longitude, point_12_4C_latitude, 'C12 (4C)', fontsize=12, ha='right')
plt.text(point_8_4C_longitude, point_8_4C_latitude, 'C8 (4C)', fontsize=12, ha='right')
plt.text(point_5_4C_longitude, point_5_4C_latitude, 'C5 (4C)', fontsize=12, ha='right')

# Set the title and labels
plt.title("Flight Route Comparison: 2C vs 4C Mercator")
plt.xlabel("Longitude")
plt.ylabel("Latitude")

# Display the legend
plt.legend()

# Show the plot
plt.tight_layout()
plt.show()

Total distance 2C Mercator: 1547.8721224618705 km
Total distance 4C Mercator: 1546.2770660005715 km

# Extract coordinates for the points in the 2C Lambert Conformal Conic (LCC) projection
point_8_2c_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[1]
point_8_2c_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[1]
point_12_2c_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[2]
point_12_2c_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[2]
point_5_2c_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[0]
point_5_2c_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[0]

# Calculate the distance between C12 and C8 in 2C Lambert Conformal Conic (LCC) coordinates
distance_12_8_2c_lambert = np.sqrt((point_12_2c_lambert_longitude - point_8_2c_lambert_longitude)**2 +
                                   (point_12_2c_lambert_latitude - point_8_2c_lambert_latitude)**2)
distance_12_8_2c_lambert_km = distance_12_8_2c_lambert / 1000  # Convert to kilometers

# Calculate the distance between C8 and C5 in 2C Lambert Conformal Conic (LCC) coordinates
distance_8_5_2c_lambert = np.sqrt((point_8_2c_lambert_longitude - point_5_2c_lambert_longitude)**2 +
                                  (point_8_2c_lambert_latitude - point_5_2c_lambert_latitude)**2)
distance_8_5_2c_lambert_km = distance_8_5_2c_lambert / 1000  # Convert to kilometers

# Total distance in 2C Lambert Conformal Conic projection
total_2c_lambert_projection_distance = distance_12_8_2c_lambert_km + distance_8_5_2c_lambert_km

# Extract coordinates for the points in the 4C Lambert Conformal Conic (LCC) projection
point_8_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[1]
point_8_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[1]
point_12_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[2]
point_12_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[2]
point_5_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[0]
point_5_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[0]

# Calculate the distance between C12 and C8 in 4C Lambert Conformal Conic (LCC) coordinates
distance_12_8_4c_lambert = np.sqrt((point_12_4c_lambert_longitude - point_8_4c_lambert_longitude)**2 +
                                   (point_12_4c_lambert_latitude - point_8_4c_lambert_latitude)**2)
distance_12_8_4c_lambert_km = distance_12_8_4c_lambert / 1000  # Convert to kilometers

# Calculate the distance between C8 and C5 in 4C Lambert Conformal Conic (LCC) coordinates
distance_8_5_4c_lambert = np.sqrt((point_8_4c_lambert_longitude - point_5_4c_lambert_longitude)**2 +
                                  (point_8_4c_lambert_latitude - point_5_4c_lambert_latitude)**2)
distance_8_5_4c_lambert_km = distance_8_5_4c_lambert / 1000  # Convert to kilometers

# Total distance in 4C Lambert Conformal Conic projection
total_4c_lambert_projection_distance = distance_12_8_4c_lambert_km + distance_8_5_4c_lambert_km

# Create a figure with larger size for better visibility
plt.figure(figsize=(12, 10))

# Plot the route from C12 to C8 for 2C Lambert Conformal Conic (LCC)
plt.plot([point_12_2c_lambert_longitude, point_8_2c_lambert_longitude],
         [point_12_2c_lambert_latitude, point_8_2c_lambert_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8 (2C LCC)", linewidth=3)

# Plot the route from C8 to C5 for 2C Lambert Conformal Conic (LCC)
plt.plot([point_8_2c_lambert_longitude, point_5_2c_lambert_longitude],
         [point_8_2c_lambert_latitude, point_5_2c_lambert_latitude],
         color='blue', marker='o', markersize=8, linestyle='--', label="C8 to C5 (2C LCC)", linewidth=3)

# Plot the route from C12 to C8 for 4C Lambert Conformal Conic (LCC)
plt.plot([point_12_4c_lambert_longitude, point_8_4c_lambert_longitude],
         [point_12_4c_lambert_latitude, point_8_4c_lambert_latitude],
         color='red', marker='o', markersize=8, label="C12 to C8 (4C LCC)", linewidth=3)

# Plot the route from C8 to C5 for 4C Lambert Conformal Conic (LCC)
plt.plot([point_8_4c_lambert_longitude, point_5_4c_lambert_longitude],
         [point_8_4c_lambert_latitude, point_5_4c_lambert_latitude],
         color='red', marker='o', markersize=8, linestyle='--', label="C8 to C5 (4C LCC)", linewidth=3)

# Annotate the points on the plot
plt.text(point_12_2c_lambert_longitude, point_12_2c_lambert_latitude, 'C12', fontsize=14, ha='right')
plt.text(point_8_2c_lambert_longitude, point_8_2c_lambert_latitude, 'C8', fontsize=14, ha='right')
plt.text(point_5_2c_lambert_longitude, point_5_2c_lambert_latitude, 'C5', fontsize=14, ha='right')

plt.text(point_12_4c_lambert_longitude, point_12_4c_lambert_latitude, 'C12', fontsize=14, ha='right')
plt.text(point_8_4c_lambert_longitude, point_8_4c_lambert_latitude, 'C8', fontsize=14, ha='right')
plt.text(point_5_4c_lambert_longitude, point_5_4c_lambert_latitude, 'C5', fontsize=14, ha='right')

# Set the title and labels for the plot
plt.title("Flight Route(Comparison of 2C and 4C Lambert Conformal Conic Projections)", fontsize=16)
plt.xlabel("Easting (Longitude in Lambert Conformal Conic)", fontsize=14)
plt.ylabel("Northing (Latitude in Lambert Conformal Conic)", fontsize=14)

# Display the legend
plt.legend(fontsize=12)

# Show the plot
plt.tight_layout()
plt.show()

# Print total distances
print(f'Total distance from C12 to C8 to C5 (2C LCC): {total_2c_lambert_projection_distance} km')
print(f'Total distance from C12 to C8 to C5 (4C LCC): {total_4c_lambert_projection_distance} km')

Total distance from C12 to C8 to C5 (2C LCC): 1460.612996865721 km
Total distance from C12 to C8 to C5 (4C LCC): 1460.5443208825818 km

# Define Mercator Projection Coordinates
point_8_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[1]
point_8_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[1]
point_12_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[2]
point_12_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[2]
point_5_mercator_latitude = flight_coordinates['mercator_y_2c'].iloc[0]
point_5_mercator_longitude = flight_coordinates['mercator_x_2c'].iloc[0]

# Define Lambert Conformal Conic Projection Coordinates
point_8_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[1]
point_8_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[1]
pt12_l_lat = flight_coordinates['2C_LCC_Northing'].iloc[2]
pt12_l_long = flight_coordinates['2C_LCC_Easting'].iloc[2]
point_5_lambert_latitude = flight_coordinates['2C_LCC_Northing'].iloc[0]
point_5_lambert_longitude = flight_coordinates['2C_LCC_Easting'].iloc[0]

# Plotting the Flight Route from C12 to C8 to C5 (Mercator and Lambert Comparison)
plt.figure(figsize=(12, 8))

# Plot the Mercator route (C12 to C8)
plt.plot([point_12_mercator_longitude, point_8_mercator_longitude],
         [point_12_mercator_latitude, point_8_mercator_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8 (Mercator)")

# Plot the Mercator route (C8 to C5)
plt.plot([point_8_mercator_longitude, point_5_mercator_longitude],
         [point_8_mercator_latitude, point_5_mercator_latitude],
         color='red', marker='o', markersize=8, label="C8 to C5 (Mercator)")

# Plot the Lambert route (C12 to C8)
plt.plot([pt12_l_long, point_8_lambert_longitude],
         [pt12_l_lat, point_8_lambert_latitude],
         color='green', marker='x', markersize=8, label="C12 to C8 (Lambert)")

# Plot the Lambert route (C8 to C5)
plt.plot([point_8_lambert_longitude, point_5_lambert_longitude],
         [point_8_lambert_latitude, point_5_lambert_latitude],
         color='orange', marker='x', markersize=8, label="C8 to C5 (Lambert)")

# Annotate points
plt.text(point_12_mercator_longitude, point_12_mercator_latitude, 'C12 (Mercator)', fontsize=12, ha='right')
plt.text(point_8_mercator_longitude, point_8_mercator_latitude, 'C8 (Mercator)', fontsize=12, ha='right')
plt.text(point_5_mercator_longitude, point_5_mercator_latitude, 'C5 (Mercator)', fontsize=12, ha='right')

plt.text(pt12_l_long, pt12_l_lat, 'C12 (Lambert)', fontsize=12, ha='left')
plt.text(point_8_lambert_longitude, point_8_lambert_latitude, 'C8 (Lambert)', fontsize=12, ha='left')
plt.text(point_5_lambert_longitude, point_5_lambert_latitude, 'C5 (Lambert)', fontsize=12, ha='left')

# Set the title and labels
plt.title("Flight Route Comparison: Mercator vs Lambert Conformal Conic(2C)", fontsize=16)
plt.xlabel("Longitude / Easting", fontsize=12)
plt.ylabel("Latitude / Northing", fontsize=12)

# Display the grid for better readability
plt.grid(True, linestyle='--', alpha=0.7)

# Show the legend
plt.legend()

# Adjust layout for better spacing
plt.tight_layout()

# Show the plot
plt.show()

import matplotlib.pyplot as plt
import numpy as np

# 1. 提取坐标：从 flight_coordinates 数据框中提取 Mercator 投影下的坐标
point_8_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[1]
point_8_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[1]
point_12_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[2]
point_12_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[2]
point_5_4c_lambert_latitude = flight_coordinates['4C_LCC_Northing'].iloc[0]
point_5_4c_lambert_longitude = flight_coordinates['4C_LCC_Easting'].iloc[0]

# 2. 计算距离：使用欧几里得距离公式计算 C12 到 C8、C8 到 C5 的距离 (4C Lambert 投影)
distance_12_8_4c_lambert = np.sqrt((point_12_4c_lambert_longitude - point_8_4c_lambert_longitude)**2 +
                                   (point_12_4c_lambert_latitude - point_8_4c_lambert_latitude)**2)
distance_12_8_4c_lambert_km = distance_12_8_4c_lambert / 1000  # 转换为千米

distance_8_5_4c_lambert = np.sqrt((point_8_4c_lambert_longitude - point_5_4c_lambert_longitude)**2 +
                                  (point_8_4c_lambert_latitude - point_5_4c_lambert_latitude)**2)
distance_8_5_4c_lambert_km = distance_8_5_4c_lambert / 1000  # 转换为千米

total_4c_lambert_projection_distance = distance_12_8_4c_lambert_km + distance_8_5_4c_lambert_km
print(f'Total distance from C12 to C8 to C5 (4C Lambert): {total_4c_lambert_projection_distance} km')

# 3. 绘制航程：使用 4C Lambert 投影下的坐标绘制航程
plt.figure(figsize=(8, 6))

# 绘制 C12 到 C8 的航程（4C Lambert）
plt.plot([point_12_4c_lambert_longitude, point_8_4c_lambert_longitude],
         [point_12_4c_lambert_latitude, point_8_4c_lambert_latitude],
         color='blue', marker='o', markersize=8, label="C12 to C8 (4C Lambert)")

# 绘制 C8 到 C5 的航程（4C Lambert）
plt.plot([point_8_4c_lambert_longitude, point_5_4c_lambert_longitude],
         [point_8_4c_lambert_latitude, point_5_4c_lambert_latitude],
         color='red', marker='o', markersize=8, label="C8 to C5 (4C Lambert)")

# 标注控制点
plt.text(point_12_4c_lambert_longitude, point_12_4c_lambert_latitude, 'C12', fontsize=12, ha='right')
plt.text(point_8_4c_lambert_longitude, point_8_4c_lambert_latitude, 'C8', fontsize=12, ha='right')
plt.text(point_5_4c_lambert_longitude, point_5_4c_lambert_latitude, 'C5', fontsize=12, ha='right')

# 设置标题和坐标轴标签
plt.title("Flight Route from C12 to C8 to C5 (4C Lambert Conformal Conic Projection)")
plt.xlabel("Easting (Longitude in 4C Lambert Conformal Conic)")
plt.ylabel("Northing (Latitude in 4C Lambert Conformal Conic)")

# 显示图例
plt.legend()

# 4. 提取 Mercator 投影下的坐标
point_8_mercator_latitude = flight_coordinates['North_Transformed'].iloc[1]
point_8_mercator_longitude = flight_coordinates['East_Transformed'].iloc[1]
point_12_mercator_latitude = flight_coordinates['North_Transformed'].iloc[2]
point_12_mercator_longitude = flight_coordinates['East_Transformed'].iloc[2]
point_5_mercator_latitude = flight_coordinates['North_Transformed'].iloc[0]
point_5_mercator_longitude = flight_coordinates['East_Transformed'].iloc[0]

# 5. 计算 Mercator 投影下的距离
distance_12_8_mercator = np.sqrt((point_12_mercator_longitude - point_8_mercator_longitude)**2 +
                                 (point_12_mercator_latitude - point_8_mercator_latitude)**2)
distance_12_8_mercator_km = distance_12_8_mercator / 1000  # 转换为千米

distance_8_5_mercator = np.sqrt((point_8_mercator_longitude - point_5_mercator_longitude)**2 +
                                (point_8_mercator_latitude - point_5_mercator_latitude)**2)
distance_8_5_mercator_km = distance_8_5_mercator / 1000  # 转换为千米

total_mercator_projection_distance = distance_12_8_mercator_km + distance_8_5_mercator_km
print(f'Total distance from C12 to C8 to C5 (Mercator): {total_mercator_projection_distance} km')

# 6. 绘制 Mercator 投影下的航程
plt.plot([point_12_mercator_longitude, point_8_mercator_longitude],
         [point_12_mercator_latitude, point_8_mercator_latitude],
         color='green', marker='o', markersize=8, label="C12 to C8 (4C Mercator)")

plt.plot([point_8_mercator_longitude, point_5_mercator_longitude],
         [point_8_mercator_latitude, point_5_mercator_latitude],
         color='orange', marker='o', markersize=8, label="C8 to C5 (4C Mercator)")

# 标注控制点
plt.text(point_12_mercator_longitude, point_12_mercator_latitude, 'C12', fontsize=12, ha='right')
plt.text(point_8_mercator_longitude, point_8_mercator_latitude, 'C8', fontsize=12, ha='right')
plt.text(point_5_mercator_longitude, point_5_mercator_latitude, 'C5', fontsize=12, ha='right')

# 设置标题和坐标轴标签
plt.title("Flight Route(Comparison between 4C Mercator Projection and 4C Lambert Projection)")
plt.xlabel("Longitude (Mercator Projection)")
plt.ylabel("Latitude (Mercator Projection)")

# 显示图例
plt.legend()

# 展示图表
plt.tight_layout()
plt.show()

Total distance from C12 to C8 to C5 (4C Lambert): 1460.5443208825818 km
Total distance from C12 to C8 to C5 (Mercator): 1546.2770660005715 km

data_geodetic = [['12 - 8 Geodetic Distance', distance_12_8_km], ['8 - 5 Geodetict Distance', distance_8_5_km], ['Total Geodetic Distance', total_distance]]
geodetic_table = pd.DataFrame(data_geodetic, columns=['Projection', 'Distance'])

geodetic_table

data = [['12 - 8 Mercator Distance', distance_12_8_mercator_2c], ['8 - 5 Mercator Distance', distance_8_5_mercator_2c], ['Total Mercator Distance', total_mercator_projection_distance_2c]]

mercator_data_table_2c = pd.DataFrame(data, columns=['Projection', 'Distance'])
mercator_data_table_2c

data = [['12 - 8 Mercator Distance', distance_12_8_mercator_km], ['8 - 5 Mercator Distance', distance_8_5_mercator_km], ['Total Mercator Distance', total_mercator_projection_distance]]

mercator_data_table = pd.DataFrame(data, columns=['Projection', 'Distance'])
mercator_data_table

data_lambert = [['12 - 8 lambert Distance', dist_12_8_l_km], ['8 - 5 lambert Distance', distance_8_5_lambert_km], ['Total lambert Distance', total_l_dist]]

lambert_data_table = pd.DataFrame(data_lambert, columns=['Projection', 'Distance'])

lambert_data_table

data_ls_lambert = [['12 - 8 LS lambert Distance', distance_12_8_4c_lambert_km], ['8 - 5 LS lambert Distance', distance_8_5_4c_lambert_km], ['Total LS lambert Distance', total_lambert_projection_4c_distance]]

ls_lambert = pd.DataFrame(data_ls_lambert, columns=['Projection', 'Distance'])

ls_lambert

speed = 185
Geodetic_travel_time = total_distance/speed
mercator_travel_time_2control_points=total_mercator_projection_distance_2c/speed
mercator_travel_time_4control_points = total_mercator_projection_distance/speed
lambert_travel_time_2control_points = total_l_dist/speed
lambert_travel_time_4control_points = total_lambert_projection_4c_distance/speed


data_travel_time = [['Geodetic', Geodetic_travel_time],['Mercator(2 control points)',mercator_travel_time_2control_points],['Mercator(4 control points)', mercator_travel_time_4control_points], ['Lambert Conformal Conic(2 control points)',lambert_travel_time_2control_points],['Lambert Conformal Conic(4 control points)',lambert_travel_time_4control_points]]

travel_time_table = pd.DataFrame(data_travel_time, columns=['name', 'time(h)'])
travel_time_table

#Bell Long Ranger Helicopter Cruise Fuel Consumption: 27 gallons per hour (gph)
#Cruise airspeed: 115 mph (185 km/h)
# Aviation Fuel Price: $1.993 per gallon
#Total Distance Flown: Total known distance flown (distance that can be calculated from Step 1)
#Step 2: Calculate Fuel Consumption
# Calculate Flight Time: We can calculate the flight time by dividing the distance flown by the cruise airspeed.
# Calculate Fuel Consumption: Fuel Consumption = Flight Time × Fuel Consumption per Hour (gph)
# Calculate Fuel Cost: Fuel Cost = Fuel Consumption × Price per gallon of fuel ($1.993)

# Assuming all the necessary variables (distances) are defined earlier

# Parameters
speed = 185  # speed in km/h
fuel_consumption_per_hour = 27  # fuel consumption in gallons per hour
fuel_price_per_gallon = 1.993  # price per gallon in USD
gallon_to_liter = 3.78541  # conversion factor from gallon to liter

# Calculate travel times
Geodetic_travel_time = total_distance / speed
mercator_travel_time_2control_points = total_mercator_projection_distance_2c / speed
mercator_travel_time_4control_points = total_mercator_projection_distance / speed
lambert_travel_time_2control_points = total_l_dist / speed
lambert_travel_time_4control_points = total_lambert_projection_4c_distance / speed

# Create the data table with travel times
data_travel_time = [
    ['Geodetic', Geodetic_travel_time],
    ['Mercator(2 control points)', mercator_travel_time_2control_points],
    ['Mercator(4 control points)', mercator_travel_time_4control_points],
    ['Lambert Conformal Conic(2 control points)', lambert_travel_time_2control_points],
    ['Lambert Conformal Conic(4 control points)', lambert_travel_time_4control_points]
]

# Create the initial travel time table
travel_time_table = pd.DataFrame(data_travel_time, columns=['name', 'time(h)'])

# Add columns for fuel consumption (gallons), fuel consumption (liters), and fuel cost
travel_time_table['fuel_consumption (gallons)'] = travel_time_table['time(h)'] * fuel_consumption_per_hour
travel_time_table['fuel_consumption (liters)'] = travel_time_table['fuel_consumption (gallons)'] * gallon_to_liter
travel_time_table['fuel_cost (USD)'] = travel_time_table['fuel_consumption (gallons)'] * fuel_price_per_gallon

# Add a new column for distance (combined distances from all projections)
travel_time_table['Total Distance (km)'] = [
    total_distance,
    total_mercator_projection_distance_2c,
    total_mercator_projection_distance,
    total_l_dist,
    total_lambert_projection_4c_distance
]

# Show the updated travel time table
travel_time_table

# Save the table to an Excel file
travel_time_table.to_excel('travel_time_table.xlsx', index=False)
processed_4c_madagascar_transformed.to_excel('processed_4c_madagascar_transformed.xlsx', index=False)

	crtl_pt_ID	Lat_deg_N	Long_deg_E	pixel_x	pixel_y	Lon_rad_E	Lat_rad_N
0	A	0.0	30.0	419	163	0.523599	0.000000
1	B	0.0	60.0	586	163	1.047198	0.000000
2	C	-30.0	60.0	586	340	1.047198	-0.523599
3	D	-30.0	30.0	419	340	0.523599	-0.523599
4	C1	NaN	NaN	525	233	NaN	NaN
5	C2	NaN	NaN	532	252	NaN	NaN
6	C3	NaN	NaN	528	253	NaN	NaN
7	C4	NaN	NaN	528	259	NaN	NaN
8	C5	NaN	NaN	517	294	NaN	NaN
9	C6	NaN	NaN	513	309	NaN	NaN
10	C7	NaN	NaN	494	301	NaN	NaN
11	C8	NaN	NaN	497	272	NaN	NaN
12	C9	NaN	NaN	499	256	NaN	NaN
13	C10	NaN	NaN	510	254	NaN	NaN
14	C11	NaN	NaN	517	242	NaN	NaN
15	C12	NaN	NaN	525	233	NaN	NaN

	crtl_pt_ID	Lat_deg_N	Long_deg_E	pixel_x	pixel_y	Lon_rad_E	Lat_rad_N	Mercator_Easting_Latitude	Mercator_Northing_Longitude
0	A	0.0	30.0	419	163	0.523599	0.000000	3.339513e+06	-7.081002e-10
1	B	0.0	60.0	586	163	1.047198	0.000000	6.679026e+06	-7.081002e-10
2	C	-30.0	60.0	586	340	1.047198	-0.523599	6.679026e+06	-3.503475e+06
3	D	-30.0	30.0	419	340	0.523599	-0.523599	3.339513e+06	-3.503475e+06
4	C1	NaN	NaN	525	233	NaN	NaN	NaN	NaN
5	C2	NaN	NaN	532	252	NaN	NaN	NaN	NaN
6	C3	NaN	NaN	528	253	NaN	NaN	NaN	NaN
7	C4	NaN	NaN	528	259	NaN	NaN	NaN	NaN
8	C5	NaN	NaN	517	294	NaN	NaN	NaN	NaN
9	C6	NaN	NaN	513	309	NaN	NaN	NaN	NaN
10	C7	NaN	NaN	494	301	NaN	NaN	NaN	NaN
11	C8	NaN	NaN	497	272	NaN	NaN	NaN	NaN
12	C9	NaN	NaN	499	256	NaN	NaN	NaN	NaN
13	C10	NaN	NaN	510	254	NaN	NaN	NaN	NaN
14	C11	NaN	NaN	517	242	NaN	NaN	NaN	NaN
15	C12	NaN	NaN	525	233	NaN	NaN	NaN	NaN

	crtl_pt_ID	Lat_deg_N	Long_deg_E	pixel_x	pixel_y	Lon_rad_E	Lat_rad_N	Mercator_Easting_Latitude	Mercator_Northing_Longitude	negative_y	mercator_x_2c	mercator_y_2c
0	A	0.0	30.0	419	163	0.523599	0.000000	3.339513e+06	-7.081002e-10	-163	3.306408e+06	-6.398864e+06
1	B	0.0	60.0	586	163	1.047198	0.000000	6.679026e+06	-7.081002e-10	-163	6.627947e+06	-6.381905e+06
2	C	-30.0	60.0	586	340	1.047198	-0.523599	6.679026e+06	-3.503475e+06	-340	6.609972e+06	-9.902338e+06
3	D	-30.0	30.0	419	340	0.523599	-0.523599	3.339513e+06	-3.503475e+06	-340	3.288434e+06	-9.919297e+06
4	C1	NaN	NaN	525	233	NaN	NaN	NaN	NaN	-233	5.407582e+06	-7.780361e+06
5	C2	NaN	NaN	532	252	NaN	NaN	NaN	NaN	-252	5.544878e+06	-8.157550e+06
6	C3	NaN	NaN	528	253	NaN	NaN	NaN	NaN	-253	5.465219e+06	-8.177846e+06
7	C4	NaN	NaN	528	259	NaN	NaN	NaN	NaN	-259	5.464610e+06	-8.297182e+06
8	C5	NaN	NaN	517	294	NaN	NaN	NaN	NaN	-294	5.242271e+06	-8.994430e+06
9	C6	NaN	NaN	513	309	NaN	NaN	NaN	NaN	-309	5.161190e+06	-9.293178e+06
10	C7	NaN	NaN	494	301	NaN	NaN	NaN	NaN	-301	4.784103e+06	-9.135992e+06
11	C8	NaN	NaN	497	272	NaN	NaN	NaN	NaN	-272	4.846716e+06	-8.558893e+06
12	C9	NaN	NaN	499	256	NaN	NaN	NaN	NaN	-256	4.888120e+06	-8.240459e+06
13	C10	NaN	NaN	510	254	NaN	NaN	NaN	NaN	-254	5.107107e+06	-8.199563e+06
14	C11	NaN	NaN	517	242	NaN	NaN	NaN	NaN	-242	5.247552e+06	-7.960179e+06
15	C12	NaN	NaN	525	233	NaN	NaN	NaN	NaN	-233	5.407582e+06	-7.780361e+06

	crtl_pt_ID	Lat_deg_N	Long_deg_E	pixel_x	pixel_y	Lon_rad_E	Lat_rad_N	Mercator_Easting_Latitude	Mercator_Northing_Longitude	negative_y	mercator_x_2c	mercator_y_2c
0	A	0.0	30.0	419	163	0.523599	0.000000	3.339513e+06	-7.081002e-10	-163	3.306408e+06	-6.398864e+06
1	B	0.0	60.0	586	163	1.047198	0.000000	6.679026e+06	-7.081002e-10	-163	6.627947e+06	-6.381905e+06
2	C	-30.0	60.0	586	340	1.047198	-0.523599	6.679026e+06	-3.503475e+06	-340	6.609972e+06	-9.902338e+06
3	D	-30.0	30.0	419	340	0.523599	-0.523599	3.339513e+06	-3.503475e+06	-340	3.288434e+06	-9.919297e+06
4	C1	NaN	NaN	525	233	NaN	NaN	NaN	NaN	-233	5.407582e+06	-7.780361e+06
5	C2	NaN	NaN	532	252	NaN	NaN	NaN	NaN	-252	5.544878e+06	-8.157550e+06
6	C3	NaN	NaN	528	253	NaN	NaN	NaN	NaN	-253	5.465219e+06	-8.177846e+06
7	C4	NaN	NaN	528	259	NaN	NaN	NaN	NaN	-259	5.464610e+06	-8.297182e+06
8	C5	NaN	NaN	517	294	NaN	NaN	NaN	NaN	-294	5.242271e+06	-8.994430e+06
9	C6	NaN	NaN	513	309	NaN	NaN	NaN	NaN	-309	5.161190e+06	-9.293178e+06
10	C7	NaN	NaN	494	301	NaN	NaN	NaN	NaN	-301	4.784103e+06	-9.135992e+06
11	C8	NaN	NaN	497	272	NaN	NaN	NaN	NaN	-272	4.846716e+06	-8.558893e+06
12	C9	NaN	NaN	499	256	NaN	NaN	NaN	NaN	-256	4.888120e+06	-8.240459e+06
13	C10	NaN	NaN	510	254	NaN	NaN	NaN	NaN	-254	5.107107e+06	-8.199563e+06
14	C11	NaN	NaN	517	242	NaN	NaN	NaN	NaN	-242	5.247552e+06	-7.960179e+06
15	C12	NaN	NaN	525	233	NaN	NaN	NaN	NaN	-233	5.407582e+06	-7.780361e+06

	crtl_pt_ID	Lat_deg_N	Long_deg_E	pixel_x	pixel_y	Lon_rad_E	Lat_rad_N	Mercator_Easting_Latitude	Mercator_Northing_Longitude	negative_y	mercator_x_2c	mercator_y_2c	East_Transformed	North_Transformed
0	A	0.0	30.0	419	163	0.523599	0.000000	3.339513e+06	-7.081002e-10	-163	3.306408e+06	-6.398864e+06	3.339513e+06	0.000000e+00
1	B	0.0	60.0	586	163	1.047198	0.000000	6.679026e+06	-7.081002e-10	-163	6.627947e+06	-6.381905e+06	6.661052e+06	1.695872e+04
2	C	-30.0	60.0	586	340	1.047198	-0.523599	6.679026e+06	-3.503475e+06	-340	6.609972e+06	-9.902338e+06	6.679026e+06	-3.503475e+06
3	D	-30.0	30.0	419	340	0.523599	-0.523599	3.339513e+06	-3.503475e+06	-340	3.288434e+06	-9.919297e+06	3.357487e+06	-3.520433e+06
4	C1	NaN	NaN	525	233	NaN	NaN	NaN	NaN	-233	5.407582e+06	-7.780361e+06	5.454904e+06	-1.381498e+06
5	C2	NaN	NaN	532	252	NaN	NaN	NaN	NaN	-252	5.544878e+06	-8.157550e+06	5.596059e+06	-1.758686e+06
6	C3	NaN	NaN	528	253	NaN	NaN	NaN	NaN	-253	5.465219e+06	-8.177846e+06	5.516603e+06	-1.778982e+06
7	C4	NaN	NaN	528	259	NaN	NaN	NaN	NaN	-259	5.464610e+06	-8.297182e+06	5.517212e+06	-1.898319e+06
8	C5	NaN	NaN	517	294	NaN	NaN	NaN	NaN	-294	5.242271e+06	-8.994430e+06	5.301982e+06	-2.595567e+06
9	C6	NaN	NaN	513	309	NaN	NaN	NaN	NaN	-309	5.161190e+06	-9.293178e+06	5.223948e+06	-2.894315e+06
10	C7	NaN	NaN	494	301	NaN	NaN	NaN	NaN	-301	4.784103e+06	-9.135992e+06	4.845236e+06	-2.737128e+06
11	C8	NaN	NaN	497	272	NaN	NaN	NaN	NaN	-272	4.846716e+06	-8.558893e+06	4.901959e+06	-2.160030e+06
12	C9	NaN	NaN	499	256	NaN	NaN	NaN	NaN	-256	4.888120e+06	-8.240459e+06	4.940113e+06	-1.841595e+06
13	C10	NaN	NaN	510	254	NaN	NaN	NaN	NaN	-254	5.107107e+06	-8.199563e+06	5.158694e+06	-1.800699e+06
14	C11	NaN	NaN	517	242	NaN	NaN	NaN	NaN	-242	5.247552e+06	-7.960179e+06	5.296702e+06	-1.561315e+06
15	C12	NaN	NaN	525	233	NaN	NaN	NaN	NaN	-233	5.407582e+06	-7.780361e+06	5.454904e+06	-1.381498e+06

	crtl_pt_ID	Lat_deg_N	Long_deg_E	pixel_x	pixel_y	Lon_rad_E	Lat_rad_N	Mercator_Easting_Latitude	Mercator_Northing_Longitude	negative_y	...	2C_LCC_Easting	2C_LCC_Northing	4C_East_Mercator_Transformed	4C_North_Mercator_Transformed	4C_geo_longitude	4C_geo_latitude	4C_geo_longtitude	4C_LCC_Easting	4C_LCC_Northing	Least Squares Fitting LCC residual
8	C5	NaN	NaN	517	294	NaN	NaN	NaN	NaN	-294	...	4.818452e+06	-3.143352e+06	5.297667e+06	-2.597039e+06	47.590771	-1.942196e+07	-0.250489	4.814277e+06	-3.143679e+06	5.749781e+06
11	C8	NaN	NaN	497	272	NaN	NaN	NaN	NaN	-272	...	4.554552e+06	-2.662252e+06	4.899877e+06	-2.159471e+06	44.017294	-1.982798e+07	-0.231681	4.552773e+06	-2.661294e+06	5.273540e+06
15	C12	NaN	NaN	525	233	NaN	NaN	NaN	NaN	-233	...	5.247102e+06	-2.069030e+06	5.456782e+06	-1.383782e+06	49.020161	-2.056871e+07	-0.258013	5.248302e+06	-2.071645e+06	5.642374e+06

	Projection	Distance
0	12 - 8 Geodetic Distance	918.748378
1	8 - 5 Geodetict Distance	552.428184
2	Total Geodetic Distance	1471.176562

	Projection	Distance
0	12 - 8 Mercator Distance	959521.693365
1	8 - 5 Mercator Distance	588350.429096
2	Total Mercator Distance	1547.872122

	Projection	Distance
0	12 - 8 Mercator Distance	954.913385
1	8 - 5 Mercator Distance	591.363681
2	Total Mercator Distance	1546.277066

	Projection	Distance
0	12 - 8 lambert Distance	911.887064
1	8 - 5 lambert Distance	548.725933
2	Total lambert Distance	1460.612997

	Projection	Distance
0	12 - 8 LS lambert Distance	911.837214
1	8 - 5 LS lambert Distance	548.707107
2	Total LS lambert Distance	1460.544321

	name	time(h)
0	Geodetic	7.952306
1	Mercator(2 control points)	8.366876
2	Mercator(4 control points)	8.358254
3	Lambert Conformal Conic(2 control points)	7.895205
4	Lambert Conformal Conic(4 control points)	7.894834

	name	time(h)	fuel_consumption (gallons)	fuel_consumption (liters)	fuel_cost (USD)	Total Distance (km)
0	Geodetic	7.952306	214.712255	812.773917	427.921524	1471.176562
1	Mercator(2 control points)	8.366876	225.905661	855.145549	450.229983	1547.872122
2	Mercator(4 control points)	8.358254	225.672869	854.264335	449.766028	1546.277066
3	Lambert Conformal Conic(2 control points)	7.895205	213.170545	806.937915	424.848897	1460.612997
4	Lambert Conformal Conic(4 control points)	7.894834	213.160523	806.899974	424.828921	1460.544321

---------------Madagascar Project------------------¶

Editor:Zhaocheng Tan¶

GUID:2964649T¶

Python Program 1---------------------------------start---------------------------------¶

Convert control point geographical coordinates to Mercator Projection grid coordinates.¶

Result A: Mercator Projection grid coordinates of control points 1-4.¶

Python Program 2-------------------------start---------------------------------¶

Estimate Digitizer to Mercator 2D similarity transformation constants using 2 control points (number 1 and number 3).¶

result B:Transformation constants estimated using 2 control points.¶

Python Program 3----------------------start-----------------------------¶

Coastline Digitiser Coordinates to Mercator coordinates using transformation constants (2 control points).¶

Result C: Coastline points in Mercator Projection grid coordinates.¶

east_transformed = a * pixel_x_co - b * pixel_y_co + t_x¶

north_transformed = b * pixel_x_co+ a * pixel_y_co + t_y¶

Python Program 4--------------------------------------------start-------------------------------------------------¶

Convert geographic coordinates to Lambert Confocal Conic (LCC) coordinates¶

Convert Mercator Projection grid coordinates to geographic coordinates.¶

Result D: Coastline points in geographic coordinates (Lat, Long).¶

Convert geographic coordinates to Lambert Confocal Conic (LCC) coordinates¶

4 control points-----------------------start--------------------------------¶

L + E = A * X¶

E = A * X - L¶

LS errors¶

x and y Least Squares Errors¶

compare to mercator error¶

Geographical distance¶

Mercator distance(2 controls point)¶

Mercator distance(4 controls point)¶

Lambert distance(2 controls point)¶

Lambert distance (4 control points)¶

Flight Route Comparison¶

Flight Route Comparison: 2C vs 4C Mercator¶

Flight Route(Comparison of 2C and 4C Lambert Conformal Conic Projections)¶

Flight Route (Comparison of 2C Mercator vs Lambert Conformal Conic)¶

Flight Route(Comparison between 4C Mercator Projection and 4C Lambert Projection)¶

the code compares the flight routes between points C12, C8, and C5 under two different map projections: 4C Lambert and 4C Mercator, calculating the distances in each projection and displaying the results on respective plots.¶

Distance:¶

Geodetic Distance¶

Mercator Distance(2 controls point)/Mercatort distance(4 controls point)¶

Lambert distance(2 controls point)/Lambert distance(4 controls point)¶

Geodetic Distance¶

Mercatort distance(2 controls point)¶

Mercatort distance(4 controls point)¶

Lambert distance(2 controls point)¶

Lambert distance(4 controls point)¶

time calculation¶

Calculation of flight time, fuel consumption and fuel costs¶

--------------------------------------end---------------------------------------¶

Student: Zhaocheng Tan 2964649T¶

Date: 15/12/2024¶